Measurement exploratory data analysis¶

In [ ]:
DEVICES = [
    'beaglebone-fan',
    'beaglebone-compressor',
    'beaglebone-pump',
    'beaglebone-refrigerator',
    'mafaulda-a',
    'mafaulda-b'
]

DEVICE = DEVICES[1]
T_WAVEFORM = 5  # (1 = MaufaulDa, 5 = others)
T_SEC = T_WAVEFORM
NFFT = 2**10   # (2**10 = MaufaulDa, 2**14 = others)
F_LIMIT = None # (3000 = MaufaulDa, None = others)
In [ ]:
import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt

from tabulate import tabulate
from IPython.display import Markdown, HTML
from tqdm.notebook import tqdm

from typing import List, Tuple
from scipy.signal import find_peaks, butter, lfilter
from tsfel.feature_extraction.features import fundamental_frequency

from zipfile import ZipFile
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda, selection, discovery, models


def beaglebone_measurement(filename: str, fs: int) -> Tuple[str, pd.DataFrame]:
    g = 9.81
    milivolts = 1800
    resolution = 2**12
    columns = ['x', 'y', 'z']
    ts = pd.read_csv(filename, delimiter='\t', index_col=False, header=None, names=columns)
        
    # Calculate amplitude in m/s^2 Beaglebone Black ADC and ADXL335 resolution (VIN 1.8V, 12bits)
    for dim in columns:
        ts[dim] = ts[dim] * (milivolts / resolution)  # ADC to mV
        ts[dim] = (ts[dim] / 180) * g                 # mV to m/s^2 (180 mV/g)
        ts[dim] -= ts[dim].mean()

    ts['t'] = ts.index * (1 / fs)
    ts.set_index('t', inplace=True)
    return (os.path.basename(filename), ts, fs, ts.columns)  # last is feature columns


def beaglebone_dataset(filenames: List[str], fs: int) -> List[Tuple[str, pd.DataFrame]]:
    dataset = []
    for filename in filenames:
        name, ts, fs, cols = beaglebone_measurement(filename, fs)
        dataset.append((name, ts))
    return dataset


def lowpass_filter(data, cutpoint, fs, order=5):
    b, a = butter(order, cutpoint, fs=fs, btype='lowpass')
    y = lfilter(b, a, data)
    return y


def mafaulda_dataset(
        place=['ax', 'ay', 'az'],
        features_path =  '../../datasets/features_data/',
        mafaulda_path='../../datasets/MAFAULDA.zip',
        rpm=2500,
        lowpass_hz=10000):

    metadata_filename = os.path.join(features_path, selection.MAFAULDA_METADATA)
    faults = {
        'normal': 'normal',
        'imbalance': 'imbalance',
        'horizontal-misalignment': 'misalignment',
        'vertical-misalignment': 'misalignment',
        'overhang-cage_fault': 'cage fault',
        'underhang-cage_fault': 'cage fault',
        'underhang-ball_fault': 'ball fault',
        'overhang-ball_fault': 'ball fault',
        'overhang-outer_race': 'outer race fault',
        'underhang-outer_race': 'outer race fault',
    }

    metadata = pd.read_csv(metadata_filename, index_col='filename')
    metadata.reset_index(inplace=True)
    metadata = models.fault_labeling(metadata, faults)
    files = pd.DataFrame()
    # Worst severity and mid rpm
    for name, group in metadata[metadata['rpm'] >= rpm].groupby(by='fault', observed=False):
        files = pd.concat([
            files,
            group[
                group['severity_level'] == group['severity_level'].max()
            ].sort_values(by='rpm', ascending=True).head(1)
        ])
    ordering = {
        'normal': 0,
        'misalignment': 1,
        'imbalance': 2,
        'cage fault': 3,
        'ball fault': 4,
        'outer race fault': 5,
    }
    source = ZipFile(mafaulda_path)
    dataset = len(files) * [0]
    for index, file in files.iterrows():
        ts = mafaulda.csv_import(source, file['filename'])
        ts = ts[place]
        ts.columns = ts.columns.str.extract(r'(\w)$')[0]
        for axis in ts.columns:
            ts[axis] = lowpass_filter(ts[axis], lowpass_hz, file['fs'])
        pos = ordering[file['fault']]
        dataset[pos] = ((file['fault'] + ' (' + file['filename'] +')', ts))

    return dataset

Load dataset

In [ ]:
if DEVICE == 'beaglebone-fan':
    Fs = 2500
    path = '../../inspections/fan/'
    files = [
        '1_still.tsv', '2_still.tsv', '3_still.tsv',
        '1_up.tsv', '2_up.tsv', '3_up.tsv',
        '1_down.tsv', '2_down.tsv', '3_down.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-compressor':
    Fs = 2500
    path = '../../inspections/datacentres/shc3/'
    files = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-pump':
    Fs = 2500
    path = '../../inspections/bvs/'
    files = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv' 
        #, 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-refrigerator':
    Fs = 2500
    path = '../../inspections/home-refrigerator/'
    files = [
        'ch1.tsv', 'ch2.tsv', 'ch3.tsv', 'ch4.tsv', 'ch5.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'mafaulda-a':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['ax', 'ay', 'az'])

elif DEVICE == 'mafaulda-b':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['bx', 'by', 'bz'])
In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

k3_1.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_2.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_3.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_4.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_1.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_2.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_3.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_4.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))
    ts.boxplot(grid=True)
    plt.show()

k3_1.tsv

x y z
count51000 51000 51000
mean 2.30926e-15 1.14551e-15 3.23924e-16
std 1.37145 2.11749 0.512991
min -2.46212 -5.23984 -1.22658
25% -1.31252 -1.59941 -0.436227
50% 0.076596 -0.377948 -0.0530239
75% 1.39386 1.32252 0.47388
max 2.68717 4.69949 1.26424
No description has been provided for this image

k3_2.tsv

x y z
count51000 51000 51000
mean -1.12739e-15 3.32283e-15 5.67347e-15
std 1.37462 1.90235 0.54732
min -2.44795 -4.07667 -1.19419
25% -1.32229 -1.4661 -0.475688
50% 0.0907731 -0.436238 0.00331546
75% 1.43198 1.07262 0.53022
max 2.48579 4.1143 1.36848
No description has been provided for this image

k3_3.tsv

x y z
count51000 51000 51000
mean 9.60069e-16 -1.59607e-15 3.9874e-16
std 1.40137 1.95904 0.55607
min -2.62483 -4.38867 -1.25534
25% -1.25967 -1.56255 -0.488938
50% 0.0336402 -0.269236 0.0140165
75% 1.4467 1.26358 0.540921
max 2.76396 4.16155 1.28338
No description has been provided for this image

k3_4.tsv

x y z
count51000 51000 51000
mean -2.11937e-15 4.49606e-15 5.62652e-15
std 0.627253 0.8348 0.22793
min -9.79898 -6.27935 -1.67302
25% -0.0273018 -0.0283453 -0.0444018
50% 0.0205986 -0.0043951 -0.0204516
75% 0.0445488 0.0435053 0.0274488
max 7.75651 9.81518 1.72791
No description has been provided for this image

k5_1.tsv

x y z
count51000 51000 51000
mean 3.58434e-15 -2.53483e-15 6.39767e-16
std 1.49955 2.32511 0.231579
min -3.22158 -5.17207 -0.665095
25% -1.09002 -2.01065 -0.162141
50% 0.251196 -0.382032 -0.0184398
75% 1.2571 2.15669 0.149212
max 2.88572 4.95886 1.25092
No description has been provided for this image

k5_2.tsv

x y z
count51000 51000 51000
mean -9.78738e-16 2.91908e-15 2.61535e-15
std 1.5101 2.32367 0.22692
min -3.3225 -4.90644 -0.651701
25% -1.09513 -2.05637 -0.148747
50% 0.270034 -0.451703 -0.0289962
75% 1.29989 2.13492 0.138655
max 2.97641 4.79339 1.19246
No description has been provided for this image

k5_3.tsv

x y z
count51000 51000 51000
mean -5.92676e-15 -7.54011e-16 1.70753e-15
std 1.49944 2.33441 0.225151
min -3.29453 -4.94197 -0.632079
25% -1.09111 -2.044 -0.153075
50% 0.298001 -0.415387 -0.00937392
75% 1.25601 2.12333 0.134327
max 2.76487 4.90156 1.21209
No description has been provided for this image

k5_4.tsv

x y z
count51000 51000 51000
mean 8.0361e-15 8.27908e-15 4.04591e-16
std 1.40058 2.10618 0.258499
min -4.15798 -5.68461 -1.33782
25% -1.02051 -1.85257 -0.140307
50% 0.0572508 0.0155409 -0.0445065
75% 1.20686 1.85971 0.123145
max 5.87715 6.93715 3.76357
No description has been provided for this image

Statistical tests

  • Normality test: Kolmogorov–Smirnov test
  • Normality visual test: Quantile-quantile plot on chosen recording
  • Stationarity test: Augmented Dickey–Fuller test
  • Stationarity visual test: Autocorrelation plot
In [ ]:
from statsmodels.tsa.stattools import adfuller
from statsmodels.api import qqplot
from scipy.stats import kstest

normality_tests = []
for name, ts in DATASET:
    for x in ts.columns:
        p_value = kstest(ts[x], 'norm').pvalue
        test = {'name': name, 'axis': x, 'p-value': p_value, 'not-normal': p_value < 0.05}
        normality_tests.append(test)

normality_tests = pd.DataFrame.from_records(normality_tests)
print(normality_tests.value_counts('not-normal'))
normality_tests.describe()
not-normal
True    24
Name: count, dtype: int64
Out[ ]:
p-value
count 24.0
mean 0.0
std 0.0
min 0.0
25% 0.0
50% 0.0
75% 0.0
max 0.0
In [ ]:
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    qqplot(ts[x], line='45', ax=ax[i], marker='.', alpha=0.5)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()
k3_1.tsv
No description has been provided for this image
In [ ]:
stationarity_tests = []
for name, ts in tqdm(DATASET):
    for x in ts.columns:
        result = adfuller(ts[x].loc[T_WAVEFORM:T_WAVEFORM+1])
        p_value = result[1]
        test = {
            'name': name,
            'axis': x,
            'statistic': result[0],
            'p-value': p_value,
            'stationary': p_value < 0.001
        }
        stationarity_tests.append(test)

stationarity_tests = pd.DataFrame.from_records(stationarity_tests)
print(stationarity_tests.value_counts('stationary'))
stationarity_tests['p-value'].describe()
  0%|          | 0/8 [00:00<?, ?it/s]
stationary
True     23
False     1
Name: count, dtype: int64
Out[ ]:
count    2.400000e+01
mean     7.193267e-04
std      3.523967e-03
min      0.000000e+00
25%      0.000000e+00
50%      1.006159e-29
75%      2.375520e-20
max      1.726384e-02
Name: p-value, dtype: float64
In [ ]:
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    ax[i].acorr(ts[x], maxlags=50)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()
k3_1.tsv
No description has been provided for this image

Time domain histogram

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].hist(figsize=(15, 3), grid=True, bins=100, layout=(1, 3), edgecolor='black', linewidth=0.5)
    plt.show()

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

No description has been provided for this image

k5_1.tsv

No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image

Time domain waveform

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

No description has been provided for this image

k5_1.tsv

No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image

Time domain waveform zoom detail

In [ ]:
for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

No description has been provided for this image

k5_1.tsv

No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image

Time domain waveform zoom - faults side by side

In [ ]:
fig, ax = plt.subplots(len(DATASET), 3, figsize=(12, 15), sharex=True)

for idx, df in enumerate(DATASET):
    name, ts = df
    columns = ts.columns
    ax[idx][1].set_title(name)
    ax[idx][0].set_ylabel('Amplitude [m/s^2]')

    for pos, axis in enumerate(columns):
        data = ts[axis].loc[T_WAVEFORM:T_WAVEFORM+0.3]
        ax[idx][pos].plot(data.index, data, linewidth=1, color='darkblue')
        ax[idx][pos].grid()
    

plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
def spectogram(x, debug=True):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))


def compare_limited_specrograms(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        #ilast = len(freqs[freqs < F_LIMIT])
        
        ax[i].plot(freqs, pxx)
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Amplitude [m/s^2]')
        ax[i].set_xlim(0, F_LIMIT)
        ax[i].set_ylim(0, 0.4)
        ax[i].set_title(name)
        i += 1


def spectrogram_energy_left_cumulative(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        
        ax[i].plot(freqs, np.cumsum(pxx) / np.sum(pxx))
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Cumulative energy [%]')
        #ax[i].set_xlim(0, 10000)
        ax[i].set_title(name)
        i += 1

Compare mafaulda faults

In [ ]:
compare_limited_specrograms(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
compare_limited_specrograms(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
compare_limited_specrograms(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image

Compare cumulative sums

In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image

Spectrogram in X axis

In [ ]:
x_spectra = get_spectrograms(DATASET, 'x')

k3_1.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_2.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_3.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_4.tsv (X axis, Fundamental = 47.7941 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_1.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_2.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_3.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_4.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in X axis

In [ ]:
show_spectrogram_detail(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

No description has been provided for this image

k3_2.tsv (X axis @ 5s)

No description has been provided for this image

k3_3.tsv (X axis @ 5s)

No description has been provided for this image

k3_4.tsv (X axis @ 5s)

No description has been provided for this image

k5_1.tsv (X axis @ 5s)

No description has been provided for this image

k5_2.tsv (X axis @ 5s)

No description has been provided for this image

k5_3.tsv (X axis @ 5s)

No description has been provided for this image

k5_4.tsv (X axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in X axis

  • MMS peak finder algorithm
In [ ]:
show_mms_peaks(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

No description has been provided for this image

k3_2.tsv (X axis @ 5s)

No description has been provided for this image

k3_3.tsv (X axis @ 5s)

No description has been provided for this image

k3_4.tsv (X axis @ 5s)

No description has been provided for this image

k5_1.tsv (X axis @ 5s)

No description has been provided for this image

k5_2.tsv (X axis @ 5s)

No description has been provided for this image

k5_3.tsv (X axis @ 5s)

No description has been provided for this image

k5_4.tsv (X axis @ 5s)

No description has been provided for this image

Harmonic series detection in X axis

In [ ]:
# show_harmonic_series(x_spectra, 'x', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.685933 1
1 285.645 0.180302 5.85
2 524.902 0.0620602 10.75
31052.25 0.0087056 21.55
4 717.773 0.0068952514.7

k3_2.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.703858 1
1285.645 0.136529 5.85
2524.902 0.0330748 10.75
3954.59 0.0057881619.55
4 14.64840.00447683 0.3

k3_3.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.656088 1
1 285.645 0.149954 5.85
2 144.043 0.0424332 2.95
31164.55 0.0068831 23.85
4 639.648 0.0053969213.1

k3_4.tsv (X axis @ 5s)

  • Fundamental frequency: 4.8828125 Hz
f y 1x
0 4.882810.0768066 1
1 153.809 0.00184415 31.5
2 297.852 0.00152948 61
3 83.0078 0.00119425 17
41015.62 0.0010988 208

k5_1.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.921647 1
1288.086 0.102319 5.9
2434.57 0.0270554 8.9
3722.656 0.0109019 14.8
4400.391 0.00575804 8.2

k5_2.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.92099 1
1288.086 0.10624 5.9
2480.957 0.0224154 9.85
3 14.64840.00968361 0.3
4627.441 0.0055030312.85

k5_3.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.936292 1
1 288.086 0.106108 5.9
2 241.699 0.0264785 4.95
3 295.41 0.0106427 6.05
41103.52 0.0058045922.6

k5_4.tsv (X axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.898246 1
1288.086 0.0986729 5.9
2480.957 0.027933 9.85
3295.41 0.00986706 6.05
4983.887 0.0044258520.15

Spectrogram in Y axis

In [ ]:
y_spectra = get_spectrograms(DATASET, 'y')

k3_1.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_2.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_3.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_4.tsv (Y axis, Fundamental = 47.3039 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_1.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_2.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_3.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_4.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in Y axis

In [ ]:
show_spectrogram_detail(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

No description has been provided for this image

k3_2.tsv (Y axis @ 5s)

No description has been provided for this image

k3_3.tsv (Y axis @ 5s)

No description has been provided for this image

k3_4.tsv (Y axis @ 5s)

No description has been provided for this image

k5_1.tsv (Y axis @ 5s)

No description has been provided for this image

k5_2.tsv (Y axis @ 5s)

No description has been provided for this image

k5_3.tsv (Y axis @ 5s)

No description has been provided for this image

k5_4.tsv (Y axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in Y axis

In [ ]:
show_mms_peaks(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

No description has been provided for this image

k3_2.tsv (Y axis @ 5s)

No description has been provided for this image

k3_3.tsv (Y axis @ 5s)

No description has been provided for this image

k3_4.tsv (Y axis @ 5s)

No description has been provided for this image

k5_1.tsv (Y axis @ 5s)

No description has been provided for this image

k5_2.tsv (Y axis @ 5s)

No description has been provided for this image

k5_3.tsv (Y axis @ 5s)

No description has been provided for this image

k5_4.tsv (Y axis @ 5s)

No description has been provided for this image

Harmonic series detection in Y axis

In [ ]:
# show_harmonic_series(y_spectra, 'y', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82811.20267 1
1478.516 0.198349 9.8
2285.645 0.103564 5.85
3573.73 0.018807411.75
4622.559 0.010690212.75

k3_2.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82811.03286 1
1 144.043 0.138609 2.95
2 285.645 0.0364966 5.85
3 541.992 0.0092398 11.1
41240.23 0.0061403225.4

k3_3.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82811.00388 1
1 478.516 0.12413 9.8
2 429.688 0.0633835 8.8
31164.55 0.0102918 23.85
41115.72 0.0077010922.85

k3_4.tsv (Y axis @ 5s)

  • Fundamental frequency: 9.765625 Hz
f y 1x
0 9.765620.0299383 1
1 21.9727 0.00412541 2.25
2 83.0078 0.00197565 8.5
3244.141 0.0015170225
4227.051 0.0013613 23.25

k5_1.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82811.38965 1
1385.742 0.147958 7.9
2529.785 0.061545910.85
3400.391 0.0299462 8.2
4 26.85550.0151213 0.55

k5_2.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82811.40826 1
1434.57 0.179166 8.9
2419.922 0.0564781 8.6
3405.273 0.0231326 8.3
4673.828 0.011014113.8

k5_3.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82811.42221 1
1 480.957 0.15385 9.85
2 578.613 0.054704711.85
3 402.832 0.0217529 8.25
41149.9 0.011927723.55

k5_4.tsv (Y axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82811.3829 1
1480.957 0.143644 9.85
2578.613 0.052124211.85
3983.887 0.016837620.15
4 14.64840.0122953 0.3

Spectrogram in Z axis

In [ ]:
z_spectra = get_spectrograms(DATASET, 'z')

k3_1.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_2.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_3.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k3_4.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_1.tsv (Z axis, Fundamental = 0.0490 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_2.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_3.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

k5_4.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in Z axis

In [ ]:
show_spectrogram_detail(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

No description has been provided for this image

k3_2.tsv (Z axis @ 5s)

No description has been provided for this image

k3_3.tsv (Z axis @ 5s)

No description has been provided for this image

k3_4.tsv (Z axis @ 5s)

No description has been provided for this image

k5_1.tsv (Z axis @ 5s)

No description has been provided for this image

k5_2.tsv (Z axis @ 5s)

No description has been provided for this image

k5_3.tsv (Z axis @ 5s)

No description has been provided for this image

k5_4.tsv (Z axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in Z axis

In [ ]:
show_mms_peaks(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

No description has been provided for this image

k3_2.tsv (Z axis @ 5s)

No description has been provided for this image

k3_3.tsv (Z axis @ 5s)

No description has been provided for this image

k3_4.tsv (Z axis @ 5s)

No description has been provided for this image

k5_1.tsv (Z axis @ 5s)

No description has been provided for this image

k5_2.tsv (Z axis @ 5s)

No description has been provided for this image

k5_3.tsv (Z axis @ 5s)

No description has been provided for this image

k5_4.tsv (Z axis @ 5s)

No description has been provided for this image

Harmonic series detection in Z axis

In [ ]:
# show_harmonic_series(z_spectra, 'z', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.280503 1
1444.336 0.0312579 9.1
2383.301 0.0217816 7.85
3507.812 0.014363710.4
4859.375 0.013331817.6

k3_2.tsv (Z axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.32505 1
1 429.688 0.0253058 8.8
2 144.043 0.0191175 2.95
3 954.59 0.016307819.55
41049.8 0.011750821.5

k3_3.tsv (Z axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.323556 1
1334.473 0.0279629 6.85
2668.945 0.021913313.7
3456.543 0.0152852 9.35
4922.852 0.010976 18.9

k3_4.tsv (Z axis @ 5s)

  • Fundamental frequency: 4.8828125 Hz
f y 1x
0 4.882810.0167918 1
1 61.0352 0.0045176712.5
2 90.332 0.0031232218.5
3324.707 0.0019177966.5
4424.805 0.0017308587

k5_1.tsv (Z axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0713598 1
1 434.57 0.0264955 8.9
2 192.871 0.0225645 3.95
3 34.17970.0192975 0.7
41010.74 0.016903720.7

k5_2.tsv (Z axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.06825791
1192.871 0.02608753.95
2336.914 0.02175146.9
3144.043 0.02002252.95
4 31.73830.01253010.65

k5_3.tsv (Z axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.8281 0.07736151
1144.043 0.02577842.95
2 7.324220.01978510.15
3288.086 0.01695075.9
4422.363 0.01405448.65

k5_4.tsv (Z axis @ 5s)

  • Fundamental frequency: 48.828125 Hz
f y 1x
0 48.82810.0907712 1
1144.043 0.0249169 2.95
2625 0.019953812.8
3964.355 0.017212519.75
4 14.64840.0139624 0.3

Histogram

In [ ]:
axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

No description has been provided for this image

k5_1.tsv

No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image
In [ ]:
axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

No description has been provided for this image

k5_1.tsv

No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image

Orbitals of all cross sections

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

No description has been provided for this image

k5_1.tsv

No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image

Orbitals of 1x harmonic frequency

In [ ]:
x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

k5_1.tsv

No description has been provided for this image
No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image
In [ ]:
t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

k3_1.tsv

No description has been provided for this image

k3_2.tsv

No description has been provided for this image

k3_3.tsv

No description has been provided for this image

k3_4.tsv

k5_1.tsv

No description has been provided for this image

k5_2.tsv

No description has been provided for this image

k5_3.tsv

No description has been provided for this image

k5_4.tsv

No description has been provided for this image